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Abstract 

We consider biochemical reaction chains and investigate how random 
external fluctuations, as characterized by variance and coefficient of varia- 
tion, propagate down the chains. We perform such a study under the assump- 
tion that the number of molecules is high enough so that the behavior of the 
concentrations of the system is well approximated by differential equations. 
We conclude that the variances and coefficients of variation of the fluxes will 
decrease as one moves down the chain and, through an example, show that 
there is no corresponding result for the variances of the chemical species. 
We also prove that the fluctuations of the fluxes as characterized by their 
time averages decrease down reaction chains. The results presented give 
insight into how biochemical reaction systems are buffered against external 
perturbations solely by their underlying graphical structure and point out the 
benefits of studying the out-of-equilibrium dynamics of systems. 

1 Introduction 



In [0Q and B2l we began a study of biochemical reaction systems subjected to 
random, external forcing. The question we considered, and continue with here, 
is the following: if we add random, external forcing to the input of a biochem- 
ical reaction system, how do those fluctuations (characterized by their variance 
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and coefficient of variation) propagate through the entire system? The broader 
aims of this paper are to gain a better understanding of how the network topol- 
ogy of biochemical reaction systems suppresses or otherwise alters the behavior 
of fluctuations in reaction systems and to point out the benefits of studying the 
out-of-equilibrium dynamics of systems. 

In [|2l we studied systems under the two simplifying assumptions that the ki- 
netics were all mass action and that each reaction involved turning precisely one 
substrate into another substrate. Therefore, in the terminology of Horn, Jackson, 
and Feinberg, each complex consisted of a single species ("ll6l lfTTl0 . Thus, we al- 
lowed reactions of the form A — > B, but not A + B — > C . These two assumptions 
caused the differential equations governing the concentrations of the species to be 
linear and so we referred to them as linear SSC (single species complex) systems. 
Considering linear SSC systems decreased some of the technical difficulties of 
the analysis while still allowing us to probe how different network structures af- 
fect the propagation of fluctuations. Under these assumptions we proved that the 
variances of fluxes decrease down reaction chains and that side reaction systems 
and feedback loops lower the variance of the flux out of reaction chains. A natural 
question is whether or not these results from [2] hold when we drop one or both of 
the simplifying assumptions. The main purpose of this paper is to demonstrate a 
biologically significant result from [|2) that does hold when we drop both the SSC 
and mass action assumptions: the variances and coefficients of variation of fluxes 
decrease as one moves down a non-reversible reaction chain. 

For an example of a reaction chain consider the following biochemical system: 
a substrate, S, enters the system at a constant rate, / > 0. This substrate then 
combines with an enzyme, E, to form ES, which is then degraded to some product 
substrate P plus the original enzyme. Finally, the product P leaves the system. If 
the concentration of the enzyme is taken to be so large as to be assumed constant 
and the reactions are non-reversible then the following graph faithfully models 
our system: 




where F x , F 2 and F 3 are functions that give the rates of the respective reactions. 
Let s, es, and p be the concentrations of S,ES, and P, respectively. Then, if the 
kinetic functions Fx, F 2 , and F 3 are functions of the reactant substrates only, the 
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differential equations governing the temporal evolution of the concentrations are 

s(t) = I - F 1 (s(t)) 

es(t) = Fi(s(t)) - F 2 (es(t)) (2) 
p(t) = F 2 (es(t))-F 3 (p(t))- 

If the functions Fi are differentiable, monotone increasing, and satisfy Fi(0) < 
I < lim^oo Fi(x), then it is easily seen that, independent of initial conditions, the 
system © will converge to the steady state (s, es,p) = (F-f F 2 ~ 1 (I), F 3 _1 (/)). 
However, if the input to the system (OQ) is allowed to fluctuate in time, then each 
concentration will fluctuate, and, hence, each flux, Fj, will also fluctuate. If the 
fluctuations are random, we can ask what the variance or coefficient of variation 
of each flux is with respect to that randomness, and how they relate. It is the goal 
of this paper to prove 

Var(Fi(s(t))) > Var(F 2 (ea(t))) > Var(F 3 (p(t))) (3) 
CV(FMt))) > CV(F 2 (es(t))) > CV(F 3 (p(t))), (4) 

where Var(-) and CV(-) represents variance and coefficient of variation, respec- 
tively. (Notice that the inequalities are strict.) 

The reaction chain given in (OQ) is an example of an SSC chain because each 
node of the network graph consists of a single substrate. In general, a reaction 
chain is any biochemical system of the following form: 

I Fx F 2 F n -i F n 

► W > > • • • ► <^n S 

where / £ M >0 is the constant input to the system, the complexes, Cj, are linear 
combinations of the substrates, and Fi : M>q — > M> are the reaction kinetics 
(where rrii is the number of distinct substrates composing complex Cj). In [|2l we 
showed that if the constant input I is replaced by the fluctuating in time random 
process I + £(t, uS), where £ (t, to) is either white noise or a mean zero, finite vari- 
ance, stationary stochastic process such that uS) > —I, and if the system © 
is a linear SSC chain, then for all i > 1, Var(Fj) > Var(F i+1 ), where the variance 
is computed according to the unique stationary measure to which the distribution 
of the species converges. In this paper, we prove that this result still holds when 
we drop the assumption that the kinetics are mass action and the assumption that 
each complex consists of a single species. The main assumption on the kinetics 
will be that they are monotone increasing in each of their dependent variables 
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(so, for example, we may consider Michaelis-Menten kinetics). In dropping the 
SSC assumption we will show that the result still holds when the complexes are 
composed of multiple species so long as each species appears in precisely one 
complex. Throughout, we will refer to systems for which complexes can by com- 
posed of multiple species yet each species appears in a single complex as MSC 
(multiple species complexes) systems. 

The goal of this paper is to prove that variances of fluxes caused by an external 
stochastic input decrease down a non-reversible reaction chain. We will show that 
if / is the average input to a reaction chain, then (once the system has reached its 
statistical equilibrium) the mean of each flux is also equal to /. Therefore, say- 
ing that the variances of the fluxes decrease down a reaction chain is equivalent 
to saying that the coefficients of variation of the fluxes decrease down a reaction 
chain. That is, equation © is equivalent to equation (HI). Because of this equiva- 
lence between the magnitudes of variances and the magnitudes of coefficients of 
variation, each result in this paper is stated in terms of variance alone and it is 
understood that each result is still valid if Var(-) is replaced with CV(-). 

Throughout, we allow external perturbations to be white noise processes or 
mean zero, finite variance, stationary stochastic processes. Considering white 
noise processes is useful because it allows one to make use of the Ito Calculus 
with which stronger results (i.e. fewer restrictions on the rate functions, Fj) can be 
proven than if you solely consider arbitrary perturbations. Also, if the input flux to 
a reaction system is perturbed by white noise, then all other fluxes are perturbed 
by mean zero, finite variance, stationary continuous processes. Therefore, one 
may construct continuous, stationary perturbations from white noise processes 
by allowing a pseudo-species to be perturbed by white noise and considering the 
output from the pseudo-species as the input to the reaction system of interest. In 
doing so, one is allowed to use the stronger white noise results as opposed to 
the stationary noise results. Thus, allowing both types of perturbations is quite 
natural. However, we point out that we do not feel the choice of external forcing 
is critical because the broader aim of this paper is to study the out-of-equilibrium 
dynamics of biochemical systems and both choices of perturbation achieve this 
aim. 

The layout of the paper is as follows. In Section |2] we consider SSC chains. 
In Section [3] we consider MSC chains. Complementing the main results are two 
important examples. In Section [2J Example 12.31 is a nonlinear chain perturbed 
by white noise for which the variance (and CV) of the species increase down 
the chain. Hence, there is no corresponding "decreasing fluctuation" result for 
the species of reaction chains. In Section [3j Example 13.81 demonstrates that the 
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assumption that each species is in precisely one complex is a necessary one. In 
both examples we use a Monte Carlo simulation to arrive at our conclusions. The 
proofs of all the results in this paper are found in Appendix lAl 

This paper is part of a larger research project in which the main biological 
goal is to understand how network topology affects how biochemical systems re- 
act to large-scale, random perturbations to their inputs. There are two distinct 
approaches we take in trying to achieve this goal. In the first, we apply random 
fluctuations to in silico representations of specific biological systems. We can 
then identify reactions, substrates, or whole subsystems that are buffered against 
the fluctuations, i.e. are homeostatic. We can then take the system apart piece by 
piece through in silico experimentation to discover the regulatory mechanisms that 
give rise to the homeostasis. In the second, we prove theorems about how random 
fluctuations propagate through relatively simple, but biologically relevant, sys- 
tems. We are interested in how these systems magnify or suppress fluctuations as 
this may give clues as to why these systems are structured as they are. In this sec- 
ond approach it is the out-of-equilibrium dynamics that is being probed in order 
to give information on the emergent properties of the system. This paper, like Q]| 
and B21, takes the second approach; for an example of the first, see |[T9ll . 

Due to the inherent randomness in the making and breaking of chemical bonds, 
biochemical reaction systems are, at their most fundamental level, modeled as 
jump Markov processes (|4||8| |f71lT2T1in~3ll ). However, if one scales up the vol- 
ume and number of molecules in a system while keeping the initial concentrations 
constant, then this intrinsic randomness becomes negligible at the scale of concen- 
trations. One is then able to faithfully model the concentrations of the substrates 
by a system of differential equations flfl4l0 . As in HI and [|2l we consider systems 
in this scaling limit. Thus, the random external forcing in this paper is on the scale 
of concentrations (and not of individual molecules) and the concentrations of the 
species are modeled by differential equations and not by discrete jump processes. 
For a more detailed comparison between the randomness in this paper and the 
inherent randomness of biochemical systems, see 0. 

2 SSC chains with random perturbations 

In this section we consider nonlinear SSC chains subjected to random perturba- 
tions. 
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2.1 The model 



A non-reversible SSC chain with a constant input is a biochemical reaction system 
with the following graphical structure: 

I F\ Fi F n -i F n 

v v v W 

> -Al > A 2 ► . . . > A n ► , 

where J > is the constant input to the system, Xj are the species (and com- 
plexes) of the system, and Fj : M> — > M>o are the reaction kinetics. If we let 
{x^ denote the concentrations of the species {Xi}, then the temporal evolution 
of x(t) is governed by the following differential equation: 

xi = I - Fi(xi) 

x 2 = Fi(xi) - F 2 {x 2 ) 

(7) 

•^n F n — \{X n — ij F n {x n ^. 

In the sequel we make the following standing assumptions on the functions Fj. 

Assumption 2.1. £ac/z Fj z'^ a reaZ valued C 1 function of[0, oo) with the following 
properties: 

a) Fi(0) = 0. 

b) For all x e R >0 , F[(x) > . 

c) lim FAx) > I. 



Note that condition|c]) guarantees that mass will not build up at any point along 
the chain so long as the input is kept at the constant value /. This assumption is 
also reasonable for systems for which the input is being perturbed by a mean zero 
random process and will be used to keep concentrations from escaping to infinity. 

We will consider two different classes of random perturbations of the input /. 
The first will be white in time while the second will be almost surely continuous in 
time. Since the kinetics, F, are defined only on the positive portion of the real line, 
it is important that any noise we consider as a perturbation to the input will never 
drive the concentrations of the species into the negative portion of the real line. 
Hence we will impose restrictions on the perturbations to ensure that the specie 
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concentrations stay non-negative at all times. Because we consider two different 
classes of perturbations of the input, we consider two different mechanisms to 
achieve this goal. 

In the case of the white in time random perturbation, we multiply the noise 
term, dB t , by a function, 9s(-), that turns the noise off if x\ approaches zero. This 
property of 6§ combined with the dynamics governing the concentration of X\, 
ensures that x\ remains non-negative for all time which, in turn, keeps all other 
concentrations non-negative for all time. The concentration of x\ is now governed 
by a stochastic differential equation, while the equations for the remaining X{ stay 
as in ©. That is, 

dx x = (I - F l (x 1 ))dt + (T0s(xi)dB(t) 
x 2 = -Fi(xi) - F 2 {x 2 ) 

(8) 

•^n F n — \\X n —. i) F n {x n ) , 

where B{t) = B(t,u>) is standard one dimensional Brownian motion, o E M>o, 
and for some small 5 > 0, 9 s (x) = 1 for all x > 5, 8 S (0) = and 9s(x) is C°° and 
monotone increasing. Since -Fi(O) = #5(0) = and / > 0, X\{t) > for all t > 
if xi(0) > 0. 

The second class of perturbations to the input considered in this paper are 
mean zero, finite variance, stationary, random processes, uS), that are continu- 
ous in time at almost every moment of time. To guarantee that the concentrations 
of the chemical substrates remain non-negative for all time we only consider per- 
turbations such that £(;£, uS) > —I for all t and u> (and so we no longer need the 
function Os(-) used in the white in time setting). In order to keep the reaction sys- 
tem away from equilibrium, we make the added restriction that for each choice 
of to, £(t,u>) is non-constant on all time intervals larger than some fixed value 
a = a(ui). We will typically write instead of oj). The almost everywhere 
continuity of £ (t) allows the possibility of isolated jumps and allows us to use a 
standard differential equation for x\ (in contrast to the Ito stochastic differential 
equation used in ©). In this case, the equations governing the behavior of the 
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concentrations are 

x 1 = I-F 1 {x 1 )+t{t) 
x 2 = -Fi(xi) - F 2 (x 2 ) 

(9) 

Xn Fn—l{Xn—l) Fn\%n) ■ 

To prove the existence of a stationary state we will further assume that the distri- 
bution of the noise's future is completely determined by its past. An example of 
such a £ (t) is a Markov processes whose future distribution depends only on its 
present value. 

2.2 Decreasing variance for SSC chains 

We are mainly interested in describing the system once it has settled into a sta- 
tistical equilibrium and any behavior that is transient in time has passed. Such 
statistical steady states are characterized by a stationary solution. A solution x*(t) 
is stationary if for any collection of times t\ < t 2 < ■ ■ ■ < t n and any s the dis- 
tribution of the vector (x*(ti + s), x*(t 2 + s), ■ ■ ■ , x*(t n + s)) is independent of 
s. When the forcing is Brownian as in ([8]), the solution is a Markov process and 
the distribution of x*(t) at any time t is an invariant measure for the associated 
Markov semigroup. An invariant measure, /x, is a measure on the state space of 
the system, R™, such that if the initial condition is chosen according to /i then 
solutions at any time t > are also distributed as fi. More precisely, if for all 
measurable A C W 1 , P(x(0) e A) = fi(A) implies that P(x(t) G A) = ji{A) for 
all t > 0, then jj is invariant to the dynamics of the system. Therefore, when the 
forcing is Brownian, a stationary solution exists. When the forcing is a stationary 
process £ (t) as in © more care must be taken to obtain a stationary solution to 
the dynamics as the solution need not be a Markov process. 

We will concern ourselves with the existence and basic properties of stationary 
solutions and invariant measures at the end of this section. First we state the 
principal result of the article and give a few numerical examples to illustrate its 
use. The following theorem is proved in Appendix IA.2I 

Theorem 2.2 (Decreasing variance down a nonlinear SSC chain). Let x*{t) be a 

stationary solution for the dynamics given in either equation ([8]) or ©. Then for 
all 1 < i < n and all t 

Var(F i (x*(t)))>Var(F i+1 (x* +1 (t))) . 
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We note that the variances of Theorem 12.21 are computed with respect to the 
choice of randomness, to, in B(t,u) or £(t,uj). That is, Theorem 12.21 gives the 
variance as an average over the choice of perturbations. In Subsection l2.3l we give 
a similar result except the variance is computed as a time average over a single 
path. In many natural settings (including those given in the examples below), the 
two notions are equivalent and, therefore, give the same intuition about fluctua- 
tions down reaction chains. 

We now give three examples where the preceding theorem holds. For the 
moment we will assume that the systems possess a stationary solution to which 
the statistics of the solutions converge as t — > oo. At the end of the section we 
will prove that the preceding assumption holds for any initial condition. 

Example 2.3 (Species variances need not decrease). Consider the following SSC 
chain with Michaelis-Menten kinetics 

10 + ae 5 { Xl )dB t F l {x x ) F 2 (x 2 ) 

— ► X\ — > X 2 — > , 

where a = 1, 5 = .001, F\(x\) = X\, F 2 (x 2 ) = We will see in Theorem 

I2.6l that this system possesses a unique invariant measure to which the statistics of 
the trajectories converge. Using Matlab to perform a Monte Carlo simulation we 
computed the means, variances, and coefficients of variation of the species and 
fluxes to be the following: 







%2 


Ftixx) 


F 2 (x 2 ) 


mean 


10 


5.18 


10 


10 


variance 


.5 


1.19 


.5 


.124 


CV 


.07 


.21 


.109 


.035 



As guaranteed by Theorem 12.21 the fluctuations of the fluxes decrease down the 
chain. However, no matter the measure we use (variance or CV), the fluctuations 
of x 2 are always greater than those of x\. Therefore, there is no counterpart to 
Theorem [2]2] pertaining to the species of an SSC reaction chain. 

To understand why the fluctuations of x 2 are higher than those of x\, consider 
the plot of F 2 (x) = Ylxj (1 + x) in Figure [2TTI The horizontal lines at x\ = 9, 10, 
and 11 represent possible fluxes into species X 2 , while the vertical lines show 
what the equilibrium value of x 2 would be corresponding to that input. While the 
perturbed system will never settle to an equilibrium, the kinetics will always be 
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Michaelis-Menton Kinetics 
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10 

8 

x*" 6 
4 
2 



"0 2 4 6 8 10 12 

X 2 

Figure 2.1: The horizontal lines represent possible inputs to species X 2 from 
species X\ and the vertical lines represent the value of x 2 that would give an equi- 
librium to the system for a given input. We therefore see that minor fluctuations 
in Fi(xi) = X\ can correspond with large fluctuations in ar- 
chiving the concentration of X 2 towards the solution of l^^lt) = X\ (t) . Therefore, 
minor fluctuations in the input to the species X 2 give rise to large fluctuations in 
x 2 . 

Example 2.4 (Continuous input and unbounded kinetics). Consider the following 
SSC chain 

10 + m F 1 {x 1 ) F 2 {x 2 ) 

— > X\ — ► X 2 — >, 

where — 10 < is a modified Ornstein-Uhlenbeck process defined in Appendix 
IB Fi(xi) = x\ and F 2 (x 2 ) = x\j{\ + x 2 ). Because both Fi(x) and F 2 (x) are 
unbounded as x — >• oo, we will see in Theorem 12.71 that the system possesses 
a unique stationary solution to which the statistics of the trajectories converge. 
Using Matlab to perform a Monte Carlo simulation we computed the means, vari- 
ances, and coefficient of variation of the fluxes to be the following: 





10 + £(t) 


Fi(xa) 


F 2 {x 2 ) 


mean 


10 


10 


10 


variance 


8 


6.8 


3.9 


CV 


0.28 


.26 


.20 



The variances and coefficients of variation of the fluxes decrease down the chain, 
as guaranteed by Theorem l2.2l 
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Example 2.5 (Continuous input and Michaelis-Menten kinetics). Consider the 
following SSC chain 



4 + £(t) F 1 (x 1 ) F 2 (x 2 ) 

— > X\ — > X 2 — >, 

where — 4 < £(t) < 4 is a modified Ornstein-Uhlenbeck process defined in Ap- 
pendix El Fi(xi) = llxi/(l + xi), and F 2 (x 2 ) = 10x 2 /(l + x 2 ). We will see 
in Theorem 12.81 that this system possesses a unique stationary solution to which 
the statistics of the trajectories converge. Using Matlab to perform a Monte Carlo 
simulation we computed the means, variances, and coefficient of variation of the 
fluxes to be the following: 









F 2 {x 2 ) 


mean 


4 


4 


4 


variance 


4.2 


3.3 


2.9 


CV 


.51 


.46 


.43 



The variances and coefficients of variation of the fluxes decrease down the chain, 
as guaranteed by Theorem 12 .21 

We now show the random dynamics given in equation © possesses a unique 
invariant measure. This invariant measure generates a stationary solution x*(t) 
when extended to paths from — oo to oo. Similarly, under some additional as- 
sumptions on £ and the Fj's, we show that the dynamics given in equation © 
possesses a unique stationary solution x*(t) (the concept of an invariant measure 
does not directly make sense for © since the dynamics are not necessarily Marko- 
vian, see Appendix IA. II) . In addition, in both settings we show that the statistics 
of the trajectories converge to those of the stationary solution x*(t) as t — > oo. In 
other words, for any x(0) and measurable A C W l , 

P{x{t) E A) -> P(x*(t) E A) = P(x*{0) E A), as t -> oo. (10) 

This means that the long time statistics of the solutions are independent of the 
initial condition and the result on the decrease of fluctuations is applicable on long 
time intervals. Of course in the setting of ([8]), x*(t) is distributed as the invariant 
measure \i so P(x*(0) E A) = fi(A). In the setting of ®, sc*(t,£) should be 
viewed as a function of the entire past of the noise. 

The next three results apply, respectively, in the three preceding examples 
to ensure the existence of an unique stationary solution to whose statistics the 
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statistics of arbitrary trajectories converge in time. The first result covers the case 
of white in time forcing while the second two apply to stationary forcing. The 
proofs of all three are contained in Appendix I A. 31 

Theorem 2.6 (Ergodicity of the SSC chain with white noise). Equation © pos- 
sesses a unique invariant measure, \i, on R n . Furthermore, the distribution of any 
solution to equation © converges to p, as t — ► oo. 

To prove the existence of a stationary solution to equation ®, we need to 
assume that the distribution of the future of the noise, £(£), is determined by its 
past (such as for Markov processes). This is made precise in Appendix IA.1I We 
also need additional assumptions on the functions Fj. We give two versions of 
these assumptions. 

Theorem 2.7. Let £ be as in Appendix IA.il Under the additional assumption 
that the rate functions Fj are unbounded as x — > oo, equation © possesses a 
unique stationary solution, x*(t). Furthermore, any solution x{t) to equation © 
converges to x*(t) as t —* oo. 

In the event that any of the F { are bounded, we need a bound on the size of 

£(*)• 

Theorem 2.8. Let £ be as in Appendix IA.il Define K = min{ lim Fi(x) — I}. 

i x— ►oo 

Under the additional assumption that —I < £(t) < M < K, for all t and some 
M < K, equation © possesses a unique stationary solution, x*(t). Furthermore, 
any solution x(t) to equation ®, converges to x*{t) as t — > oo. 

In the white in time setting, the system is in fact ergodic and hence by Birkhoff 's 
ergodic theorem we know that for almost every realization the time average of any 
statistic converges to the value of the statistic in the invariant measure. Combining 
this with the strong mixing properties of such a system we have that 

lim - / (F i (x i (s))-I) 2 ds = Var(F i (x*(t))) 

t^oo t Jq 

for almost every realization of the Brownian forcing and every initial condition 
Xq. For this to hold in the setting of ©, we need to assume in addition that the 
stationary measure on £ is ergodic. Even without this assumption, the next section 
shows that one can say something in general. This underlines the fact that the 
decrease of variance is really a pathwise phenomenon due to the dynamics. 
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2.3 Pathwise perturbations 



The variance described in the previous sub-section is computed with respect to the 
probability measure of the perturbations. More precisely, if to is the realization of 
the perturbation then Var(F i (x{t 1 to))) = K UJ (F(x i (t, lu)) — I) 2 , i.e. is an average 
over the realizations of the noise. Another natural way to characterize how pertur- 
bations propagate down chains is to consider the time averages of paths. Consider 
again the dynamics given by equation ©, except now the only assumptions on 
are pathwise assumptions: 

If 1 If 1 
i.) lim - / £,(s)ds = and n.)limsup- / £(s) 2 ds < oo. (11) 

t^OO t Jq t ~*O0 t Jq 

That is, we now assume that the time average for converges to zero and that 
the time average of the square is bounded above. The following theorem states 
that the pathwise variances of the fluxes do not increase down reaction chains and 
is proved in Appendix IA.4I 

Theorem 2.9. Consider equation © where satisfies (fTT|) . Then for all i > 1, 
the following hold: 

i r* 

1. lim - / Fi(xi(s))ds = I. 

t^OO t Jq 

2. liminf (- [ £(s) 2 ds-- [ (FAx^s)) - I) 2 ds) > 0. 

t^OO \t Jo t Jq J 

3. liminf Q j\Fi( Xi {a)) - Ifds - ~ j\F i+1 (x i+1 (s)) - I) 2 ds^J > 0. 



3 MSC chains with random perturbations 

We now consider MSC chains with random perturbations. We will again allow 
perturbations that are white in time or that are stationary, mean zero, finite vari- 
ance and continuous for almost every t and that satisfy the conditions of Appendix 
IA.1I Consider a reaction chain, ©, where each complex, Ci, consists of nii unique 
species and no species is contained in more than one complex. Thus, if x(t) is the 
vector representing the species concentrations at time t, then x(t) E IR mi+ '" +m ™. 
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Let Xf represent the jth species in complex i and be the multiplicity of species 
Xf in Complex i. For example, if the reaction chain is 

— ► Xl + 2Xl — ► 3X 2 X — 

then v n = 1, f 12 = 2, and v 21 = 3. 

If Fi represents the reaction rate from complex Cj to complex C i+ \ we have 
that, 

F i {x(t)) = F i (x 1 i ,...,x? i ):W ni ^R. 
We assume each Ft satisfies As sumption 13.11 which is analogous to Assumption 

Assumption 3.1. Fi is a real valued C 1 function of [0, oo) m » with the following 
properties: 

a) If x\ = Ofor any 1 < j ' < rrii, then Fi(x) = 0. 

b) If I < j <miandx E R>£, then -jjF^x) > 0. 

c) 3M > such that if xl > M for all species in the ith complex, then Fi(x) > I. 

As in the SSC case, if we want to add a random perturbation to the input flux 
of the system, we must only consider perturbations that will never drive concen- 
trations into the negative portion of the real line. We handle this issue in a similar 
manner as in the SSC case: 

If the perturbation is white in time, we multiply the perturbation by a function 
which will go to zero if the concentration of one of the species in the first complex 
goes to zero. Therefore, let 9s(x\,..., x™ 1 ) : R mi — > M> satisfy the following 
three properties for some small 5 > 0. 

1. 0g(x) — 1 when each x{ > 5. 

2. 9s(x) = if x{ = for any 1 < j < mi. 

3. 9s is C°° and is monotone increasing in each of the variables x\, . . . , x™ 1 • 

If we add a white noise perturbation multiplied by 9$ (x) to the input of the system, 
then the dynamics are now governed by the stochastic differential equation 

dx(t) = f(x(t))dt + a9 s (x)dB(t)u, (12) 
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where u = [v n , v X2 , . . . , V\ mi , 0, . . . , 0] T , a G M >0 , -B(t) is standard one dimen- 
sional Brownian motion, and f[(x) = Vu(I — Fi(x)) for each 1 < i < m x , 
/2 = Uai(-Pi(a:) - ^(x)), for each 1 < i < m 2 , . . . , jft = v„i(F n _i(a;) - F n (a;)) 
for each 1 < i < m n . 

If the perturbation is a mean zero, finite variance, stationary process, to), 
that for each uj is continuous for almost all t and that satisfies the conditions of 
Appendix |A.1[ then in order to keep the concentrations non-negative, we again 
assume that to) > —I for all t and u. In this case the dynamics of the system 
are governed by the differential equation 

= /(*(*)) +£(*)«, (13) 

where / and u are as above. 

The following four theorems are analogous to those in the SSC case and their 
proofs can be found in Appendix I A. 51 

Theorem 3.2 (Decreasing variance down a nonlinear MSC chain). Let x* (t) be a 

stationary solution for the dynamics given by either equation (TT2l) or (fT3l) . Then 
for all i > 1 and t > 

Var(Fi(x*(t))) > War {F i+1 {x* (t))) . 



Theorem 3.3 (Ergodicity of the MSC chain with white noise). Equation (PT21) 
possesses a unique invariant measure, fi, on IR n . Furthermore, the distribution of 
any solution to equation © converges to \i as t — > oo. 

Theorem 3.4. Let £ be as in Appendix IA.il Under the additional assumption 
that the rate functions F^ are unbounded as x — ■> oo, equation (fT3l) possesses a 
unique stationary solution, x*{t). Furthermore, any solution x(t) to equation (fT3l) 
converges to x*(t) as t — > oo. 

Theorem 3.5. Let £ be as in Appendix \A.l\ Define K = min{ lim F{(x) — I}. 

i x—*oo 

Under the additional assumption that —I < £(t) < M < K, for all t and some 
M < K, equation (fT3l) possesses a unique stationary solution, x*(t). Further- 
more, any solution x(t) to equation (TT3l converges to x*(t) as t — > oo. 

Example 3.6 (Sum of Two Species with Mass Action Kinetics). Consider the 
following MSC chain with mass action kinetics 

/(*) F, F 2 F 3 

— ^ Y — > X 1 + X 2 — . X 3 + X 4 — 
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where I(t) = 10 + 26 s (y)dB t (with 5 = .001), F 1 {y) = y, F 2 (x 1 ,x 2 ) = x x x 2 , 
F 3 (x 3 , x 4 ) = X3X4. Using Matlab to perform a Monte Carlo simulation we com- 
puted the means, variances, and coefficient of variation of the fluxes to be the 
following: 





Fi(y) 


^2(^1,^2) 


F 3 (x 3 ,Xi) 


mean 


10 


10 


10 


variance 


2 


1.73 


1.62 


CV 


0.14 


.131 


.127 



We note that the variances and coefficients of variation of the fluxes decrease down 
the chain, as guaranteed by Theorem 13 .2[ 

Example 3.7 (Sum of Two Species with Michaelis-Menten Kinetics). Consider 
the following MSC chain with Michaelis-Menten Kinetics 

I(t) F 1 F 2 F 3 

— > Y — ► Xi + X 2 — > X 3 + X4 — >, 

where I(t) and Fi are as in Example [331 and F 2 (xi,x 2 ) = 14xiX 2 /[(l + £i)(l + 
x 2 )}, F 3 (xz, X4) = 14x3X4/ [(1 + x 3 )(l + X4)}. Using Matlab to perform a Monte 
Carlo simulation we computed the means, variances, and coefficient of variation 
of the fluxes to be the following: 





Fi(y) 


^2(^1,^2) 


^3(^3,3:4) 


mean 


10 


10 


10 


variance 


2 


.72 


.49 


CV 


.14 


.085 


.07 



As guaranteed by Theorem 13.21 the variances and coefficients of variation of the 
fluxes decrease down the chain. 

Example 3.8 (Species can not be in more than one complex). Consider the fol- 
lowing MSC chain subjected to white noise perturbations for which the species 
Xi appears in two complexes (and so this system is not covered by Theorem 13 .21) 

I(t) Fx F 2 F 3 

— ► X\ + X 2 — > X 3 — > X\ + Xi — 

where I(t) = 10+9 s (xi,x 2 )dB(t) (with 5 = .001), F 1 (x u x 2 ) = 2x^ 2 , F 2 (x 3 ) = 
x 3 , and F 3 (xi, x 4 ) = 5xiX 4 . We performed a Monte Carlo simulation using Mat- 
lab to compute: 
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Fi(x 1 ,x 2 ) 


F 2 (x 3 ) 


F 3 (xi,x i ) 


mean 


10 


10 


10 


variance 


4.16 


.45 


1.71 


CV 


.204 


.067 


.131 



Note that Var(F 3 (x)) > Var(F 2 (x)) and CV{F 3 (x)) > CV(F 2 {x)). Therefore, 
the assumption in Theorem 13.21 that each species is in precisely one complex is 
necessary. 

4 Discussion 

We have proven under a variety of different contexts that if the input to a non- 
reversible biochemical reaction chain is perturbed by a random process, then the 
variances and coefficients of variation of the fluxes will decrease as one moves 
down the chain. The assumptions made on the different choices of perturbations 
and on the properties of the rate functions were varied, and explicitly spelled out. 
Further, much care was taken to state precisely what is meant by "fluctuations de- 
crease down reaction chains." Due to this (necessary) mathematical detail, how- 
ever, it is easy for the over-riding point of the paper to be lost: considering the 
out-of-equilibrium dynamics of a biochemical system can be an important tool 
for understanding the dynamical properties of that system. 

A comparison of the results of this paper to metabolic control analysis (MCA) 
( lfl2llfT0ll ) sheds light on the importance of considering out-of-equilibrium dy- 
namics. The control coefficient for the flux out of a reaction chain, F, in terms of 
the input, /, is 

r r_dF I 

where the values are computed at equilibrium. However, at equilibrium, F = I. 
Therefore, independent of the choice of reaction kinetics or the length of the chain, 
Cf = 1 . This implies that changes in the output of a chain correspond directly 
with changes to the input. However, by studying the out-of-equilibrium dynamics, 
we have shown in this paper that the fluctuations in a reaction chain will actually 
decrease as one moves down a reaction chain and changes to the output of a 
chain do not correspond directly with changes to the input. The differing results 
are biologically significant since it is tempting to speculate that this decrease in 
fluctuations (and, hence, increase in stability) is one reason long reaction chains 
may be evolutionarily advantageous in cellular systems. 
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While all of the technical details of the proofs have been relegated to the ap- 
pendices, we would like to point out that to prove the main results of this paper 
(with the exception of Theorem |2.9l) two things must be shown: 1) the existence of 
a unique solution whose statistics are stationary for the dynamics and 2) variances 
of a stationary solution decrease down reaction chains. Typically, a stationary 
solution can be proved to exist so long as the perturbations to the system do not 
drive any solutions to infinity. Stability properties of the non-perturbed system 
can then be used to show uniqueness of the stationary solution. The fact that the 
variances decrease down reaction chains follows from standard inequalities and 
the use of Lyapunov type functions. Intuitively, however, the variances decrease 
down reaction chains because the dynamics are always forcing the output flux 
from a complex towards the input flux. That is, the dynamics are constantly mov- 
ing the system towards a shifting equilibrium. There is, however, a natural time 
delay in its ability to do so. Therefore the output will always be lagging behind 
the input, which leads to the decrease in variance. 

There is still much work to be done in studying biochemical reaction systems 
subjected to external perturbations. A natural extension of this work and that 
of [0 would be to attempt to analyze reaction systems with more complicated 
geometries and more complicated kinetics (like product inhibition). The main 
technical issues encountered in such a study would be: 1) the extremely weak 
stability of many such systems ([5][6]) would make proving the existence of a 
stationary solution difficult and 2) it will be difficult to isolate the variances of 
particular fluxes or specie concentrations within a complicated system. While 
both of these problems are formidable in a theoretical study such as in this paper, 
they become trivial in an in silico study ( |fT9l ). 
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A Precise definition of the noise f and the proofs 



A.l Assumptions on the noise £ (t) needed for existence 

In addition to the standing assumptions that £ (t) is stationary with mean zero and 
finite variance, to prove the existence of a stationary solution for ©, we need 
to assume that the distribution of £ (t) is determined entirely by the past of £ on 
any interval of time (— oo, s] with s < t. Intuitively, we mean that given the 
value of £(s) for s £ (— oo, t) the distribution on [t,T] is uniquely determined 
for any T > t. This must be done in a way such that if one first adds a segment 
[t, T — r] and then [T — r, T], the resulting distribution on [t, T] is the same as 
if one had added the segment [t, T] in one step. For a discussion of some of the 
issues involved if one does not make such an assumption, see [|9l and subsequent 
works by the author. 

Let C ae ((— oo, 0], R) denote the space of almost every where continuous func- 
tions, /, endowed with the norm sup \f(s)\e~ a<1 ^ for some a > 0. Let V t be a 
Markov semigroup on C ae ((— oo, 0], R) which is Feller, has an invariant measure 
Ai, and such that for A'f-almost every 70 £ C ae ((— 00, 0], R), Vt{jo, ■ ) is con- 
centrated on elements ^ t £ C ae ((—oo, 0], R) with 7t(s) = 7o(s + t) for s < — i. 
If 7t 5 t > 0, is a realization of the Markov chain generated with V t with 70 dis- 
tributed as AL we define £(s) = 7 t (s — t) for s < t. This is well defined since 
our assumptions on Vt make 7 s (r — s) = 7 t (r — i) for r < min(s, t). 

The dynamics of such a P t can be understood as follows. Given an initial 
history from —00 to 0, one adds on a segment of length t, resulting in a trajectory 
from — 00 to t. After shifting this trajectory back by —t, one again obtains a 
trajectory from —00 to 0. The distribution of this new trajectory from —00 to 
is given by Vt- The conditions above simply insure that the trajectory from —00 
to —t coincide with the initial trajectory from —00 to zero. If is a Markov 
process, then it can be constructed as above and hence is an example of the type 
of noise we allow. 

A.2 Proof of principal result on variances 

Proof of Theorem 1272] We consider the dynamics given by equation ©. The 
proof when the dynamics is given by equation © is identical. 



19 



Defining ^(t) = Fi (#*(£)) — /, the equations governing x\ and x 2 are 

dx\ = (I- F 1 {x* 1 ))dt + a6 s {xl)dB{t) (14) 
x* = F 1 (x* 1 ) - F 2 (x* 2 ) =1- F 2 (x* 2 ) + 6- (15) 

We claim that for any t, EF^xKt)) = EF 2 (x^(t)) = I. Integrating CH), taking 
expected values, using that the distribution of x*(t) is stationary and noting that 
a L 9s(xl(s))dB s is an L 2 — martingale gives 

Exl(t) = EarJ(O) + (/ - EFi(x*(t))) t. 

By the stationarity of the system, Ex$(t) = Ea:*(0), so E {Fi(x*(t)) - I) = as 
claimed. A similar argument which uses EFi (a;* (t)) = I (and, hence, E£i(£) = 0) 
shows that EF 2 (x 2 (t)) = I. Therefore, in order to show that War(F 2 (x 2 (t))) < 
Var(F!(x*(t))), we need E (F 2 (x*(t)) - I) 2 < E^(t) 2 . 
Let Gi(ar) = 2 fi(F 2 (y) - I)dy. Then, 

jG x {xi{t)) = G'^ximm 

= 2(F 2 (x* 2 (t)) -W- F2te(t)) + (16) 
= -2(F 2 (x* 2 (t)) - I) 2 + 2(F 2 (x* 2 (t)) - I)Ut). 



Pick t > arbitrarily. Integrating (1161) up to time t and taking expected values 
gives 



EG 1 (x 2 (t))-EG 1 (x 2 (0)) = -2 / E(F 2 (x 2 (s)) -ifds 

Jo 

+ 2 f t E[(F 2 (x 2 (8))-I)^(a)]d 8 
Jo 



(17) 



Using that x 2 (t) is stationary, differentiation of (fT71) together with the inequality 
2ab < a 2 + b 2 gives 

= -2E(F 2 (x* 2 (t)) - I) 2 + 2E [(F 2 (x*(t)) - J)£ 
<-E(F 2 (^(i))-/) 2 + Ee 1 (t) 



We claim, however, that the inequality in (1181) is strict. To see why, we suppose, 
in order to find a contradiction, that 2E [(F 2 (x* 2 (t)) - /)&(*)] = E(F 2 (x 2 (t)) - 
J) 2 + E£i(F) 2 . Then F 2 (x* 2 (t)) - I = fi(t) with probability one. However, this 
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implies F 2 (x 2 (t)) = Fi(x*(T)) with probability one. Because t was arbitrary, 
we conclude that with probability one F 2 (x 2 (t)) = Fx(x*(t)) for all t in some 
countably dense subset of R. However, by the continuity of the functions involved, 
this implies that with probability one F 2 (x 2 (t)) = for all t G R. Thus, 

i;|(t) = for all time and x 2 (t) is a constant. But, F 2 (x* 2 (t)) = Fi(x*(t)) and 
so x\(t) is also a constant. However, for any t > 0, P{x*(s) = const : s G 
[0, £]} = 7^ 1. Thus, the inequality in (fT8l) is strict, which was the desired result 
of the Theorem. Therefore, the result is shown for the first step in the chain. To 
complete the proof, one simply repeats the argument down the chain. □ 

A.3 Existence and uniqueness of stationary solutions and in- 
variant measures 

The proofs of Theorems 12.61 [2771 and !2.8l have the same overall structure. We use 
the assumptions over the dynamics to obtain a uniform in time bound on some 
statistic of the concentration vector which can be used to prove that a sequence 
of time averages is tight. By extracting a convergent sub-sequence we can prove 
the existence of at least one invariant measure for the white in time setting. For 
the stationary forcing, we must work on the space to trajectories stretching back 
to negative infinity and prove the existence of a stationary measure on that space. 
We then prove that the invariant measure or stationary solution is unique and that 
the statistics of any solution converge to it under the dynamics of the system. 

To prove the needed tightness for Theorem 12 .61 we make use of the following 
Lyapunov function: 



where the V^'s are positive numbers yet to be determined and the Xj are defined 
as the solution to Fj(xj) = I (that is, they are the equilibrium values of the 
unperturbed problem). As an example, for a chain with n = 2 we have 



We begin by proving a fact that, while technical, is the crux of the proof of Theo- 
rem [2T61 




(19) 



V\ o V 2 2 

V{x) = —{xi- Xi) + — \{x x - Xt) + (x 2 - x 2 )] . 
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Lemma A.l. Let A be the generator of the SDEs ©. Then there are positive 
numbers Vi, V 2 , . . . , V n and positive numbers c, k such that ifV(x) is defined by 
(fT9l) then AV(x) < c - k\x\. 



Proof. For all k < n: 

dV d 

dx k dx k ^ 2 

8=1 



E ( X 3 ~ 

.i=i 



i=k 



E ( x i ~ 
.i=i 



(20) 



Let F = I. Using equation (|20l) , it can be shown that 

n ot r n / n 

fc=l j=l \ i=j 

Therefore, 

1 <9 2 " r)V 

AV{x) = r %(x i r-^V(x) + E ^r(Fk-i - F k ) 



k=l 



dx k 



\Ue 5 {x x fj2v\ + j2 
V i=i / j=i 

^(x^E^ +E 
i=i / i=i 



where the last equality is a definition. We now choose the V/s recursively. Let 
V n = 1. Because lim F n (x) > /, s n is bounded by: 

X^rOO 

Sn\XJ {x n Xfi)\I l'n\Xri) ) *^ C n k n X n , 

where c n and fc n are some positive constants. Then s n _i is given by 

Sn-l(x) = (X„_i - X n _i) - F n _i(x n _ x )) + (I - i 71 „(x n ))) . 

^(xn) > 0, so if x n _i > then 

S„-i(x) < (x n -i - X n -l) {V n -i{I ~ F„_i(x n _i)) + I) . 

We may therefore choose V n -\ to be large enough so that there are positive con- 
stants c n _i and k n -\ such that 
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Continuing up the chain, we consider Sj for j < n. When Xj > Xj we have 

W \V j (I-F j (x j ))+ J2 V ' {l 



s j\ x ) ~ \ x j x j) 

i=j+i 

< fa - Xj) (vtf - Ffa)) + 1 £ V t 

\ i=j+l 

Since Vj + i . . . , V n have already been defined, we may choose Vj so large that there 
are positive constants c, and kj such that 

s j ^ c j ~ kj x j- 



x\ 



Setting 

- n n 

c= -a 2 J2Vi + J2c h 

i=l i=l 

we now have that for some k 

n 

AV(x) < c — kiXi < c — k 
i=i 

which was the desired result. □ 

Proof of Theorem \Z6\ The proof has two parts. First, we will use Lemma lATTI 
and Prohorov's theorem ([3], pg. 59) to show that there exists a measure which is 
invariant to the stochastic flow generated by equation ([8]). We will then prove that 
this invariant measure is unique and that all distributions converge to it under the 
flow of the SDE®. 



Parti 

Let V(x) be defined by (fT9l where Vi, V2, . . . , V n are given by Lemma lATl 
Then, if k, c > are the constants given in the conclusion of Lemma lATl 

dV(x) = AV(x)dt + dM(t) < (c - k\x\) dt + dM(t), 

where M(t) is some L 2 - martingale. Integrating gives 

V(x(t)) < V(x(0)) + ct-k [ \x(s)\ds + M(t) - M(0), 

Jo 
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where x(0) is some fixed value. Rearranging terms, taking expected values and 
using the fact that EM(t) = EM (0) then yields, 

- / E|rr(s)|ds <H v v . 

t Jq K rvt 

Thus for any R > Chebychev's inequality gives 

1 /* c 1 V(x(0)) 1 



(A P {Ws) |>fl }(is <-- + ^-, (21) 

where, again, the initial condition x(0) is fixed. The right side of equation (I2TT) 
converges to zero uniformly in t > 1 as R — > oo. Therefore, the sequence of 
measures on M. defined by 



1 



u n {A) = — / P{a;(s) G 

where t n — > oo as n — > oo and A C M", is tight ([3J, pg. 59). By Prohorov's 
theorem, z/ n is relatively compact and so there exists a subsequence v nk and a 
measure /u, such that v nk — > //, where the convergence is weak convergence. 
Thus, for all A C R 



1 /" tn fc 

/i(A) = lim — / P{x(s) G A}ds. 



For A C IR n , let (j) T {A) = {x(T) : x(0) G A}. To show that // is invariant to 
the flow of equation ([8]) we need to demonstrate that for all T > and Acl, 
/40y X )(A) = where /i(0y 1 )(A) = //(x : (p T (x) G A). Note that, by 

definition, fj,(x : 4>t{x) E A) = lim^^ ^- / * nfc P{x(s + T) G A}ds. Using a 
change of variable, we then make the following computation for any T > and 
A C R: 

//(^(A) = lim — / P{x(s + T) G A}cis 



*n fc +T 



lim — / P{x(s) G A}ds + lim — / P{x(s) G A}ds 

s— >oo t nfc Jq k^oo t nk Jt„ k 

1 f T 

-lim — / P{x(s) G A}ds 

I P{x(s) G A}ds - / P{x(s) G A}ds 
Ju., Jo 



fi(A) + lim - 
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However, 





rtn 


lim — 









P{x{s) e A}ds - / P{x(s) e A}ds 



2T 

< lim — = 0, 



k— >oo t 



and so (x{4> T l )(A) = ^(A). Thus, /i is invariant under the stochastic flow gener- 
ated by ®. 



Part II 

The proof that the invariant measure is unique is not completely straightfor- 
ward. The noise enters only one species, hence the diffusion is not uniformly 
elliptic (so arguments such as in [fTSll do not suffice). The proof we now sketch 
follows a rather standard line of argument. We refer the reader to lfT6l [|20l ifTTll for 
the missing details. The proof has three elements. First, one shows that the gen- 
erator of the diffusion satisfies Hormander's "sum of squares" theorem and hence 
is hypoelliptic. This ensures that the Markov transition density pt(x, y) is smooth 
in x and y and hence is a Strong Feller process. This gives the local smoothing 
needed to ensure that the invariant measure found above is unique. The structure 
of © and the fact that the F[ do not vanish ensures that the span of the needed 
Lie brackets is of full dimension. Hence Hormander's theorem holds. 

Secondly, we need to provide the global information which ensures open set 
irreducibility (the fact that processes starting from different initial points have 
nonzero probability of entering a small neighborhood of each other). The Lya- 
punov function given by (fT9l) shows that the processes return to a bounded ball 
B about the origin eventually. Since there is a globally attracting fix point, if the 
noise is small for long enough all of the points of B will enter an arbitrarily small 
neighborhood of the fixed point. 

Finally, with the above facts in hand, the uniqueness and convergence result 
follows from standard arguments (see HI [|20l 011 II22T ). 

□ 

Proof of Theorem \277\ As in the proof of Theorem l2.6[ the proof is split into two 
parts. In the first we prove the existence of a stationary solution x*(t) for the 
dynamics ©. In the second we show that the if x(t) and y(t) are solutions driven 
by the same noise, then y{t) — > x{t) pathwise. Hence, we conclude there can only 
be one stationary solution since any two would converge to each other over time. 
Parti 

Unlike the previous example, the process x(t) alone is not a Markov process. 
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However, if we include the entire history of £ then the system does become Marko- 
vian. More precisely, from the assumptions in section |AT| we know that £(s) = 
j t (s — t) for s < t where jt is a Feller Markov process on C ae ((— 00, 0], R) with 
semi-group Vt and with 70 distributed as the invariant measure M.. Then the pair 
(x(t), 7*) is a Markov process on the expanded state space R n x C ae ((— 00, 0], R). 
Let Vt denote the Markov transition semi-group of this system and ir x and 7r 7 the 
projection onto the x and 7 coordinate respectively. Since we start j t from an 
invariant measure for its dynamics, we know that the statistics of 74 are constant 
in time equal to M. for all t and hence is tight. 

Let be an arbitrary initial condition for x(t) Defining the measure 

Qt{-) = \J^ J ^(^ (0) '7, -)M{di)ds 

we need only show that Qt^x 1 i s tight to conclude that Q t is tight since Qt^ 1 = 
M. is independent for t. We will do this coordinate by coordinate. Consider the 
equation governing x 1 (t): 



x 1 (t) = I-F 1 (x 1 (t))+£(t). (22) 



Integrating (|22|) gives 



Xl (t) =x(0) + It- [ F 1 {x 1 {s))ds+ [ £{s)ds (23) 
Jo Jo 

<x(0) + lt- / F 1 (xi(s))l { | a:(s) | >i i } ds+ / £(s)ds 
Jo Jo 

<x(0) + It-FtiR) [ l {lMs)l>R} ds + [ £{s)ds. 

Jo Jo 

Taking expected values and rearranging terms gives 

i/V {Ms ,|>^ s <^ + |if. ,24) 



Note that rearranging equation (1231) and taking expected values gives us the addi- 
tional bound 

~ [ t EF 1 (x 1 (s))ds<I+^^. (25) 
t Jo t 

Continuing down the chain we consider X2(t): 

x 2 {t)=F l {x 1 {t))-F 2 (x 2 {t)). 
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Integrating gives 

x 2 {t) = x 2 {0)+ [ F^x^ds- [ F 2 {x 2 {s))ds (26) 

Jo Jo 

<x 2 (0)+ / F 1 (x 1 (s))ds-F 2 (R) I l {lx2(s)l>R} ds 

Jo Jo 

Rearranging terms as before, taking expected values and using equation (1251) gives 

1 /"* Dri / k , „> , . / Exi(O) Ex 2 (0) 



( , o p{|*MI > *}«. < F2(R) . m(R) . m{R) 

Further, rearranging equation (|26l) and using equation (|25l) gives 

1 /"* „ , , NN , Exi(O) Ex 2 (0) 

' EF 2 (x 2 (s))ds < / + — ^+ - 



We may continue down the chain in a similar manner and conclude that there are 
positive constants c\, c 2 , . . . , c n such that for alH > 1 

1 f* 

QtT^x (iV : SU P \Vi\ > R }) = T / p { su P \ X i( S )\ > R } ds ( 27 ) 

< 



J2]J Q P{\x l (s)\>R}ds<J2 



Fi(R)' 



Because each Fj is monotone and unbounded the right side of inequality (1271) 
above converges to zero uniformly in t as R — » oo. Therefore, just as in the proof 
of Theorem 12 .61 we may invoke Prohorov's Theorem to guarantee the existence 
of a measure \i on W 1 x C ae ((— cxd, 0], K.) that is invariant to the dynamics induced 
by Vf By Kolmogormov's extension theorem we can use this measure to define a 
measure on pairs of noise £ and solution trajectories x starting at — oo and contin- 
uing to oo. The projection of this measure onto the solution coordinate produces 
a stationary solution for the x(t) dynamics. One should really view this stationary 
solution x* along with its noise trajectory £ which was constructed along with it. 
Part II 

Let x*(t) = x*(t,£) be the stationary solution and matching noise trajectory 
found above. Let y(t) be the solution starting from an arbitrary initial condition 
y(0) using the same noise £(£). 

Consider x*(t) and yi(t). If x*(0) = j/i(0), then x*(t) = yi(t) for all time by 
uniqueness of solutions. Suppose that x*(0) > yi(0) (if x\{Q) < j/i(0) there is a 
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symmetric argument). Then, x\{t) > yi(t) for all time. Differentiating x\{t) 
2/i(t) gives 

j t {x\(t) - yi (t)) = -(Fx(x;(t)) - F^it))) 



Defining 



we have that 



1 11 } 1 " x\(t)- yi (t) 



H l{ t) = jUgC^j^fl), (28 ) 



<(*) -yi(*) = (xl(0)- yi (0))e- J « Hlis)ds . 

By the uniform bound given in equation (|27l) . we know that both and yi(t) 
spend a positive fraction of time in a compact set on which H\(t) > di > 
for some di > (since -f/i is an approximation to the derivative of F{). Thus, 
— 2/1 (t) — > 0, as t — > oo. 
We next consider x^t) and j/ 2 (t). Suppose that ^(O) < 2/2(0). Let r 2 be the 
first time x* 2 (t) = 2/ 2 (i). Then, up until time r 2 , 

|(y 2 (t) - x*(t)) = -(^(xj) - F 1 (y 1 )) - (F 2 (y 2 (t)) - F 2 (x* 2 (t)))- 
But, x\ > yi so (F 1 (x* 1 ) - F^)) > and 

j t (y 2 (t) - x* 2 (t)) < -my 2 {t)) - F 2 (x* 2 (t))) 

y 2 (t) -x* 2 (t) 

Defining H 2 (t) as we did H\{€) we may conclude from the above that up until 
time r 2 



y-2 



(t) - x* 2 (t) < (2,3(0) - x* 2 (0))e-ti H *U*: 



Therefore, if r 2 = 00, then, as in the previous case, we may use the above equation 
and the bound (1271 ) to conclude that \x 2 (t) —y 2 (t) | — > 0, which is the desired result. 
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If r 2 is finite, then for all time after r 2 , x 2 (t) > y 2 {i). To see this note that if 

^(O = 2/2(i), then 

±( x ;(t)-y 2 (t))=F 1 (xl)-F 1 (y 1 )>0, 

where the inequality follows since x\ > y\. Thus, we consider times past r 2 and 
redefine our initial condition to be the values x(r 2 ) and y(r 2 ). 
We note: 

j t (x*(t) - y 2 (t)) = - (F 2 (x* 2 ) - F 2 (y 2 )) + F^) - F^y,) 

^ F 2 (x* 2 ) - F 2 (y 2 ) (29) 
= ~{x 2 - y 2 ) — + F 1 (x 1 ) - F 1 {y 1 ). 

To gain control over the term Fi(x^) — -Fi(yi) we use the equations governing x\ 
and y ± : 

x\(t) - Vl (t) = xl(0) - yi (0) + / F^xKs)) - F 1 (y 1 (s))ds. (30) 

Jo 

Rearranging and using that x\ > y x , we have 

F x [x\{s)) - F 1 (y 1 (s))ds < xl(0) - Vl (0). (31) 



J 

Jo 



Thus, if 77! (t) = J* F 1 (x* 1 (s))-F 1 ( Vl (s) ) ds, we have that F 1 (x\ (t) ) - F 1 ( Vl (t) ) 
r)[(t), and that for all t, ^(t) < a^(0) - yi(0). Therefore, 

f t (x*(t) - y 2 (t)) = -{x* 2 - y 2 )H 2 {t) + r,[(t), 
and integrating by parts gives 

x*(t)-y 2 (t) = (x*(0)-y 2 (0)e-ti H ^ dr + f V [(s)e~ ^ H ^ dr ds 

Jo 

= 04(0) -y 2 (0)e-^ Ha ^ ds + Vi{t) - Vi(0'K f « H2(r)dr 
77i(s) e -£ Ha(r)dr if 2 (s)ds. 



The last two terms are negative and, as before, the exponential terms go to zero as 
t — > 00. So, lim^oo \x 2 (t) — y 2 {t)\ < r]i(t) < 2^(0) — yi(0). However, we can 
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re-scale time (do the above analysis on the interval [t/2, t] instead of [0, t\) to con- 
clude that lim^oo \x* 2 {t)-y 2 {t)\ < lim^^i^) < lim^^ {x\{t/2) - j/i(t/2)) = 
0. 

Now we continue down the chain in a similar manner and consider x 3 and 
y 3 . Without loss of generality we may assume x 2 (t) > y 2 (t) for all time. If 
%3(t) < 2/3 W f° r ai l time, we do the same argument as above to conclude that 
\x3(t) — 2/3 ) | — > 0, as t — > oo. Thus, we assume x$(t) > ys(t) for all time. 
The argument is the same as that above, except we now have to get control over 
F 2 {x* 2 {s))-F 2 {y 2 {s)). We have 

xl(t)-y 2 (t)=x* 2 {0)-y 2 (0)+ [ - F^s))) ds 

Jo 

~ [ (F 2 (x* 2 (s))-F 2 (y 2 (s)))ds. 
Jo 

Rearranging gives 

/ (F 2 (x*(s)) - F 2 (y 2 (s))) ds = x* 2 (0) - y 2 (0) + Vl (t). 
Jo 

We may define the above integral to be rj 2 {t) and perform the same analysis as 
before. In this way we continue down the chain and conclude that lim^oo \x(t) — 
y(t)\ = 0, which was the desired result. 

□ 

Proof of Theorem^ Let N{ = F{ 1 (I + M + e), where e < K-M. If x x {t) > 
N±, then by the mono tonicity of F\ we have 

xi(t) = I-F( Xl (t))+Z(t) 

<I-F(N*)+Z(t) 
= -M-e + £(t) 
< -e. 

Therefore, independent of initial conditions, lim^oo X\(t) < N{. However, e was 
arbitrary, so lim^oo X\(t) < F^ X (I + M). Continuing in this manner down the 
chain shows lim^oo Xi(t) < Ff 1 ^ + M), for each i. Thus, for large t, there 
exists L > such that E|x(t)| < L. By Chebychev's inequality we then have 

±J*P{\x 1 (s)\>R}ds<^, 
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which converges to zero uniformly in t as R — > oo. As in the proof of Theorem 
12.71 we need to consider the Markov process on the extended state space R n x 
C ae ((— oo, 0], R). As before, we obtain tightness by using the above estimates 
on the marginal of this measure in the x(t) variable since the £(t) variable is 
stationary and hence already tight. We may again use Prohorov's Theorem to 
guarantee the existence of an invariant measure. The proof of uniqueness is the 
same is in the proof of Theorem 12 .71 □ 



A.4 Proof of Theorem 2.9 



Proof of Theorem |Z91 We begin by showing that lim^oo xi(t)/t = 0. Consider 
the dynamics governing Xi where satisfies (ITTb 

x 1 = I-F 1 (x 1 ) + £(t). (32) 

Letxi = F _1 (J). Then 

j t {x x {t) -x x ) = i- F 1 ( Xl (t)) + m = - I -—^^-M t ) - so + m 

Setting H(t) = ^——t^ > (which is well defined since F\ is assumed differ- 
entiable and is positive by the monotonicity of F{) and using Duhamel's formula 
gives us 

( Xl (t) - x x ) = (xi(0) - x 1 )e-^ H ^ da + f e-^ H{r)dr i(s)ds. 

Jo 

Integrating by parts gives 

pt 

(x 1 (t)-x 1 ) = (x 1 (Q)-x 1 )e-ti H ( s ) ds + / ^ s )ds 

Jo 

+ J e-tt H{r)dr H(s) f| \(r)dr^ ds. 
By the positivity of H(t) and property (fTTI) . we then have 

l im = i im I f ' e -SsH^ dr H{s) ( f\(r)dr) ds. 

t^oo t t^oo t Jo \J J 
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Let e > 0. There exists an S > such that s > S implies |^ J* £(r)dr\ < e/2. 
There exists a T = T(s) > such that t > T implies sup s<5 1 1 J 8 ^{r)dr | < e/2. 
Therefore, if t > max{S, T} we have that 



e -^ H{r)dr H{s) (J £{r)dr^J ds 

< f e -£ H V*H{s) 
Jo 

<e / e-^ H{r)dr H(s)ds 
Jo 



£(r)dr 



vl{s<s} + t{s>s})ds 



< e. 

Thus, ]im t _ +DO a; 1 (t)/t = 0. 

Integrating equation (|32l) . dividing by £ and taking the limit as t — > oo now 
gives us 

1 /■* 

lim - / F{x(s))ds = I, 



which proves part 1 of Theorem 12 .91 for x\. 

Let G(x) = 2 J 1 (Fi(y) — I) dy, which, non-coincidentally, is the same func- 
tion used in the proof of Theorem 12 .2[ We have 

jGMt)) = -2(F 1 (x 1 (t)) - I) 2 + 2(F 1 (x 1 (t)) - I)t(t). 
Integrating and using the inequality ah < (l/2)a 2 + (l/2)6 2 gives 

G(xi(t))<G(xi(0))- [\F 1 (x 1 (s))-I) 2 ds+ [ £(s) 2 ds. 

Jo Jo 

Therefore, part 2 of Theorem |2.9l will be shownfora;i if liminf^oo G(xi(t))/t > 
0. We have 

liminf -G(xi(t)) = 21iminf - / (F x (y) - I) (l{ y > Sl } + U y <xi}) dy 



> 2 lim inf 

t^oo t 



> -21 lim ^ 

t^oo t 



(Fi(y) - I) l [y <x 1 }dy 
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so part 2 is shown for x\ . Note that parts 1 and 2 of Theorem 12.91 show that 
F\{x\{t)) — I satisfy condition (fTTT) . Therefore, to prove parts 1, 2, and 3 for 
all Xi, one simply continues down the chain considering Fi(xi(t)) as the external 
perturbation of Xi + ±. 



A.5 Proofs of Section |3] 

The proofs of Theorems 13.21 [331 13 .41 and |3.5l can be handled simultaneously. 

Proof of Theorems \3. 21 \3.3\ \3.4\ and \3.5\ The key to each proof is the recognition 
that the species in each complex satisfy constant multiples of the same (stochastic) 
differential equations. Mathematically this means these species can by grouped 
and treated as a single substrate with a redefined kinetics. This reduces us to the 
case previously studied. More explicitly, there are constants and dyk such that 
xl(t) = dijkXi(t) + Cijk for all t. Thus, the species of each complex can be solved 
for from knowledge of just one species from that complex. Further, a monotone 
increase in one translates to a monotone increase in the others. Choosing one 
species, yi, from each complex, we may redefine the Fj's (and 9$ in the white 
noise case) appropriately so that the vector function y{t) satisfies either © or 
© with the F/s satisfying Assumption 12.11 Therefore, applying the theorems of 
Section I2T21 completes the proof. □ 



B Processes Used in Examples [23] and 123 



In Example |2.4| £ (t) is described as a modified Ornstein-Uhlenbeck process such 
that —10 < More precisely, is governed by the following dynamics: 



dm 



-£(t)dt + 4dB(t) if f (t) > -10 



-i{t)dt 



if £(*) < -10 



This dynamics ensures that if £(0) > —10 then £(£) > —10 for all t. 

In Example 12.51 £ (t) is built from the Ornstein-Uhlenbeck equation d£ (t) = 
-£{t)dt + 3dB(t), with the added condition that if f (t) = -4 or £(t) = 4, then 
d£(t) = —£(t)dt. More precisely, is governed by the following dynamics: 



-£(t)dt 

d£(t) = { -£{t)dt + ZdB{t) 

-mdt 



if i{t)dt > 4 

if - 4 < £(t) < 4 

if £{t)dt < -4 
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Figures with Captions 

Figure I2TT1 

Michaelis-Menton Kinetics 

12 
10 
8 

>T 6 
4 
2 



2 4 6 8 10 12 

X 2 



Caption: 

The horizontal lines represent possible inputs to species X 2 from species X x 
and the vertical lines represent the value of x 2 that would give an equilibrium to the 
system for a given input. We therefore see that minor fluctuations in Fi(xi) = X\ 
can correspond with large fluctuations in x 2 . 
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